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Abstract 

We present a Hamiltonian approach for the wellknown Eigen model of the 
Darwin selection dynamics. Hamiltonization is carried out by means of the 
embedding of the population variable space, describing behavior of the sys- 
tem, into the space of doubled dimension by introducing additional dynamic 
variables. Besides the study of the formalism, we try to interpret its basic 
elements (phase space, Hamiltonian, geometry of solutions) in terms of the 
theoretical biology. A geometric treatment is given for the considered system 
dynamics in terms of the geodesic flows in the Euclidean space where the 
population variables serve as curvilinear coordinates. 

The evolution of the distribution function is found for arbitrary distributed 
initial values of the population variables. 
PACS: 87.10.+e 
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I. INTRODUCTION 



Methods of Hamiltonian and Lagrangian analysis are extensively applied 
during last decades in various fields of theoretical and mathematical physics. 
These methods, however, have not been widely used until now in mathemat- 
ical theory of biological and ecological systems which foundations originated 
from the works of the 10-30-es by A. Lotka, V. Volterra, J. Haldane and R. 
Fisher. 

The basic obstacle, by our view, is that the mathematical models of bi- 
ological subjects, being open systems in their essence, are not supposed to 
be Hamiltonian ones from the outset. On the other hand, broad potential- 
ities may appear for biological systems to analyze some aspects (of both 
principle and technical character) of their dynamics, if the system admits 
Hamiltonization in a sense. 

In the present paper we suggest a Hamiltonian form for the wellknown 
Eigen model of the Darwin selection dynamics. To construct the Hamiltonian 
formalism, we extent the space of population variables introducing additional 
degrees of freedom. As a result the behavior of the system is described in 
terms of the space with doubled dimension if compared to the original one. 
Besides the study of the formalism, we try to interpret its basic elements 
(phase space, Hamiltonian, geometry of solutions) in terms of the theoretical 
biology. A geometric treatment is given for the considered system dynamics 
in terms of the geodesic flows in the Euclidean space where the population 
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variables serve as curvilinear coordinates. 

The notion of a Darwin system (DS) as a formal object in theoretical 
biology had been introduced by M. Eigen to deduce laws of living system 
evolution from the principles of theoretical physics and chemistry. 

By definition, the DS is an open system, which is constituted of coupling 
units of different species (genotypes) self-copying with a small number of er- 
rors (convariantly self-reduplicating according to Timofeev-Resovskii nomen- 
clature [jlj]). The coupling units utilize a substance and free energy of an 
external nutrient supply |2j . 

The Darwin selection phenomenon occurs in such systems under the fol- 
lowing necessary conditions: the stability of system organization (i.e. the 
total quantity of units of all species should conserve) or the stability of feed 
component influx. Both of these constraints imply the flow through the 
system. A simple experimental model of the DS (chemostat, turbidostat, 
etc.) can be a system where viruses and bacteria are defined as convariantly 
self-reduplicating units (see, for example, ||). 

The following system of differential equations proposed by M. Eigen to 
describe the DS evolution: 

N 

±i = x i (A l Q l - Di) + J2 w ji x j ~ (!) 

here k = 1,...,N; N is a species quantity of self-reduplicating units 
(genotypes) in the system; X{ is a specific quantity (concentration) of the i-th 
species; A{ (Di) is a specific reproduction (death) velocity of the i-th species; 
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Qi is a parameter of the reduplication quality of the i-th species (takes the 
values from to 1); Wjj is a specific velocity of the error reduplication of 
the j—th species into the i— th one; ^ is a dilution parameter which usually 
equals to Dxi. Here D is a dilution flow rate. 

Further, as it has been shown by M. Eigen ||, the system ([I]) can be put 
in a simpler form 



using obvious relations between Ai, Qi, wji and by introducing new pop- 
ulation variables (yi), being related to the so called " quasispecies" which 
are described by certain combinations of original variables (xj). Here 
i, j,k = 1, N; N is a (constant) quasispecies quantity in the system; y^ 
has a meaning of a specific quantity (concentration) of the i-th quasispecies; 
S = (Si, . . . , Sf) are the concentrations of the external nutrient supply com- 
ponents; fii(S) is a generalized Maltuzian parameter which has a meaning of 
the reproduction specific velocity of the i-th quasispecies and depends on S. 
The fii(S) is an algebraic combination of all A^, Dk,wjk in which Ai makes 
a major contribution. The fii(S),D have dimension of inverse time. The 
"quasispecies" introduced by M. Eigen as new population variables have a 
clear biological sense of clones, systems of organisms with a prevalence of 
the determined genotype in the system and small admixtures of some other 
genotypes connected with the leading one by the mutation and recombination 
transitions @. 




(2) 
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There are two commonly recognized types of DS under the constraints 
imposed above on the system ([I]). Recall the definitions of these types. 

Definition 1.1 A Darwin system is referred to as DS with stable organiza- 
tion (DSSO) if E^=i Hi and S are constant. A Darwin system with D = const 
is referred to as DS with a stable flow (DSSF). 

It is known that the representation of the system (|2J) in terms of relative vari- 
ables having a meaning of shares (probabilities) of the quantities || trans- 
forms (0) into the classical form of the Fisher system of equations which 
describes the Darwin selection dynamics in the panmictic populations under 
the special requirement for the fitness function [|5j. The system (g) is also 
valid for description of interspecies competition for the grow-determining 
substrate when there is no migration. In view of significance of the systems 
(IID^O) m biology and ecology they are investigated in a large number of 
works since the paper j5j has appeared. So we touch upon some of these 
publications relevant for the present paper. 

An exact solution of (g) is known only for /ij = const j6j], a general form 
of approximate solution is found by Jones j7j] for the \±i changing in time. 
A general review of the properties of system Q2j) in the Fisher modification 
is given in |8j, |9j. Stability of the stationary solutions of (g) is analyzed 
in detail by Pykh fpTOfl - Feistel and Ebeling have studied behavior of the 
Eigen-Fisher system on the adaptive Right's landscape and shown that in 
a number of cases the system can be represented in the canonical gradient 



form , [jnj . Isida []T2j] investigated non-equilibrium thermodynamics of the 



systems described by the Eigen equations. As is known, there are two dif- 
ferent approaches to describe Darwin, Fisher and reducible to them Volterra 
systems by means of extremal principles. The first (teleological) approach 
assumes to extremalize of a functional having a certain biological sense as 
a finite aim of the evolution. For example, it might be maximum of an av- 
erage population fitness [j5|], or maximum of an average productivity [2j, or 
maximum of a reproductive potential j8j] etc. The second approach assumes 
to achieve the extremum of some function of dynamic variables in each mo- 
ment of time. Thus, in reference [O] it is shown that the requirement for 



Kulback's measure to have a maximum value leads to the relations which are 
the solutions of (g) under the special conditions imposed on the system in 
the approximation of the DSSO. As a next example we refer to the work |9|] 
where variational form of the Eigen-Fisher dynamics in a special case is built 
up by introduction of a metrics of a Riemann space. Despite the fact that 
the dynamics of the Eigen-Fisher type systems are managed to represent in 
the form of extremal principles in some special cases, the problem stays far 
from complete solution as it has been explained in j8|]. 

The Hamiltonian form of the DS dynamic equations, being proposed in 
the present work, leads to the conventional variational principle, which is 
habitual in mechanics. 
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II. DS WITH A STABLE ORGANIZATION 

In accordance with def.fTT], Eqs. (@) assume the following form for the 
DSSO case: 

/ N . N \ 

yi = y% yn - E vm / E %j > ^ = const - ( 3 ) 

The dimensionless variables Pi = Vi/ E^Li 2/j have a meaning of shares (prob- 
abilities) of the quantities. Taking into account that EjLi 2/j = const, we 
derive from (g) for ^ the equations: 

N 

Pi=Pi{tH-Y,Pjf J 'j)- ( 4 ) 

.7 = 1 

If j9tv is eliminated solving the constraint Y,f=iPj = 1, one can rewrite (f|) 
as follows: 

Pi =Pi(Oi - E ( 5 ) 

Here and up to the end of this section we assume: i, j, k = 1, . . . , N — 1, a* = 
Hi — Hn, is a Maltuzian parameter of iV-th quasispecies taken as a gauge 
one, for example, by the criterium of minimal value of jii. Next, to simplify 
(0) we introduce the following variables: 

N-1 

Zi = log]pi/(l- E^)]- (6) 

3 = 1 

Accounting the condition 1 — Y,j=i Pj = Pn, the new variables Z{ might 
be called as informational ones, inasmuch as they reflect, in a spirit of the 
Shannon's theory [[T|j , a quantity of information for every independent degree 



of freedom co-ordinated to the eliminated N-th variable. 



Consider variables (§) as coordinates of a manifold A which naturally can 
be named as an information one. A state of the system in A is characterized 
by a point z = (zi, z% . . . , Zjv-i)- The system dynamics z(t) is described by a 
velocity vector z(£). The system (|5|) in the variables z takes the simple form, 

k = a u (7) 

and has the general solution 

Zi{t) = Z iQ + (Lit. (8) 

Hence, evolution of the DSSO in the "information" variables (0) is a uniform 
rectilinear motion. Taking "mass" to be a unit, we can conceive that the ai 
is the i-th component of the momentum. In accordance with Hamiltonian 
formalism, we will interpret (zi) and (a^) as the position and momentum 
variables, respectively, in a 2(N — 1)— dimensional phase space M. of the 
system. Introducing Hamiltonian of the system 

I N-l 
1 3=1 

we can write down the Hamiltonian form of the system (0) as follows: 

« = dk/da, = a,, (9) 
di = —dh/dzi = 0. 

In the framework of the considered Hamiltonian approach we can state the 
following variational principle for the DSSO dynamics. Lagrange function 
of the system (0), / = EjLl 1 ^^— h, after the momentum variables (a^) 
eliminating, takes the form: 



1 

i j=1 



(10) 



The minimal action principle is: 



A = f l(z,z){r)dr, 5 A = 0. 



(11) 



The variational principle (|IU|) , ( |TT|) leads to the following treatment of the DS 
dynamics: 

a DSSO evolutes in the information manifold A in such a way which minimal- 
izes the squared norm of the respective velocity of variation of information 
in the system at the considered time interval. 

Canonical form of the dynamic equations (|5|) is found from the canonical 
transformation (zi,di) — > (pi,Ti), where are additional variables canonically 
conjugate to pi. This transformation is: 



Pi = e*/(1 + A), Asit 1 ^ 



(12) 



r J = (l + A)[fl i e-H$ 1 o,i 
The inverse transformation: 

e*=pi/(l-n), a^Ef^ 1 ^, 
a* = p«(ri - iu), it; = E^Ll 1 r^pj. 
Eqs. (|9|),([T^),([T3|) lead to the original system (|5|) to be supplemented by the 
equations for n: 



(13) 



Pi = Pi[Pi( r i - - EjLi 1 p]( r J - W)], 
h = r i Sjtl 1 Pj{ r j - w ) - Pi( r i ~ m ) 2 - 



(14) 



Using (g), (|I2f),(|I3l) we obtain the solution of ( IHQ in the form 



/,n exp[p 0l {r 0l - w )t] f 

Pi{t)=P0i JfZi > ( 15 ) 

1 - Q + I] Pok exp[p ok (r ok - w )t] 
k=l 



n(t) = {r 0i exp(-poi(r j - w Q )t) + w [l - exp(-p i(roi - w )t)]} x 

(16) 

{l - Q + E£iVofc exp(p fe(roA; - w )t)} . 
Here j^o/c = Pk{t)\t=o, rot = fk(t)\t=o are the initial conditions for the system 
0, w = Ek=iPokrok, = Ek=iPok- Expressions r (t,p,r), Po(t,p,r) are 



obtained from (|I5|), fllfil) by the substitution: po« ^ Pi? fo* ^ fi, t —> —t. 
The system ( \P$j ) is the Hamiltonian one in terms of the canonical coordinates 



iPh r i) {Pi are the position variables and T{ are the momentum ones) with the 
Hamiltonian 

1 N-l 

h = - Y,p){r j -wf. (17) 

1 3=1 

The respective dynamics can be represented in terms of the geodesic flows 
of the (N — 1)— dimensional Euclidean space where pi serves as curvilinear 
coordinates and z% are the Cartesian ones. To impart tensor nature to the 
notations, let Z{ — > z\ pi — > p l leaving subscripts at The Hamilton 

function (JTTJ) can be rewritten as 



1 N-l 

^=^E 9 %3 {p)nr r (18) 

2 i=i 

Here g lJ (p) are the contravariant components of the metric tensor in the 
coordinates (pi): 
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9 ij (p) = E 1 = W 2 S ij ~ W) V + 2?V) 2 ] + (19) 

& — 1 

N-l 

where R 2 = E (j 9 *) 2 - The covariant components are: 

k=l 

. N N ~ 1 dz k dz k Sij 1 1 iV- 1 

P«(P) = E IT-IT- = + 1771 m + 



fc ti <9p*<V (p*) 2 ^(l-ft) ^'(1-ft) (i-^) 2 ' 



det(^(p)) = (i - ^nJLiV) 2 , 



(20) 



Y,k=\ 9 %k 9kj = ^j, {bij,^ are the Kronecker deltas). The Lagrange func- 
tion, 

N-i i N-l 

I = E Afc - h = - E 9ki(p)p k p l , 

k=l z k,l=l 

defines the Euler-Lagrange equations, 

d/dt(dl/dp>) - dl/dp> = 0, (21) 

which take the form of geodesic equations: 

i2 j n-i . 

-Et+Y. n,{v)v k v l = o, 

ai k,l=l 

1 N-l 

where T J kl (p) = - E 9 JS {dg ks / dp 1 + dg is /dp k - dg M /dp s ) are the Cristoffel 

2 s =i 

symbols. 

Let us point out that the coordinates pi are nonorthogonal ones due to 
nondiagonality of the metrics in these coordinates. 

The DS trajectories pi(t) are usually mapped as curved lines in the coor- 
dinates pi supposing them to be orthogonal Cartesian ones (see, for example, 
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), whereas this motion can be represented as straight lines in the Cartesian 
coordinates zi. Twisting of the trajectory with respect to (p^ is the result of 
the curvilinearity of these coordinates. The existence of the coordinates (2^), 
where the DSSO dynamics is represented by a uniform rectilinear motion, 
illustrates the absence of a coupling between quasispecies and any external 
factors in the course of the Darwin selection process. From this standpoint, 
the DSSO can serve as a biology-theoretic analogue of an inertial reference 
system in the Newtonian mechanics. 

Let us note that the a\ give rise to an Abelian algebra of integrals of the 
system ([P|) which is integrable by Liouville with respect to the conventional 
Poisson brackets {/, g} = J2j=i(df/dpjdg/drj— dg/dpjdf /drj). The system 

dpi/da = da q /dr l = p q (5 qt - pi), 

{22) 

dr t /da = -ddq/dpi = -5 iq {r q - w) + p q r t 

defines a flow generated by a q in the phase space M. which leaves Eqs. (|T3j) 
invariant, a is a group parameter. Integration gives: 

p q (a) = p q ($)e a V(a)-\ 
Pi(a) =Pi(0)V(a)- 1 , i^q; 

(23) 

r q (a) = e- a V(a) 2 [r q (0) + W(e a - l)V(a)- 1 ]] 
n(a) = ri(0)V(a), i ^ q. 

Here V(a) = 1 + (e a - l)p q {0), W = ^(0)^(0), ^(0), n(0) are initial 
conditions for (|22|) . Let us note, that the variables pi in (|23| ) transform 



independently of r« so that one can directly use for the original system 
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In conclusion of this section we consider DSSO dynamics with less num- 
ber of constraints then in the def.|TTT|. Let us retire only the requirement of 
constant of total quantity of all species in the system and remove the restric- 
tions on the nutrient supply concentration and on some other parameters of 
the medium, for example, on a temperature. This leads for \i{ to be certain 
functions of time. Since fii are not involved in the transformation from the 
original variables to the "information" ones, Hamiltonian form of the DS 
dynamics can be thereby constructed in the considered case too. 

Let fi be the first derivative of hi = fi. Since ai serve as momentum 
variables in the DS dynamics, then the derivatives of fi can be considered 
as components of a force. It is easy to show that the Hamiltonian of the 
considered system takes the form: 

N I 

h(z,a) = Y,(^-z l f l ). (24) 

Respectively, Hamiltonian form of Eqs. (0) is: 

Zi = dh/dai = ai, 

(25) 

h t = -dh/dz l = fi. 

The Hamiltonian ( gjj) plays a role of an "energy" in the Darwin dynamics 
if functions fi are constant in time. The second term in ( gj| ) has a sense of 
a "potential energy". The DS dynamics described by the system ( |25| ) is a 



motion in the field of the "force" fi. The additional "potential" term leads 
to the variational principle flTUp,(|TTf) with the modified Lagrangian: 
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III. DS WITH A STABLE FLOW 



To provide Eqs. ([I]) to be a complete system in the case of DSSF (in the 
sense of def. pl]) , one has to define explicitly the functions /i^(S) and to extend 
the system introducing equations for S. Here and further i,j,k= 1, . . . , N. 

The wellknown Monod function |T5{] is commonly recognized for DSs. It 



depends on one independent variable S which has a meaning of concentration 
of the substrate being in physiological minimum (growth limiting substrate). 
That is the concentration vector S is reduced to one component S. The 
dynamic system in this case takes the form: 

Vi = Vi{^i{S) - D), 

S = D(So-S)-Ef =1 y^j, (27) 

Hi = m l S/(K l + S), 
where rrii is a maximal reproduction specific velocity, K{ is a saturation 
parameter of z-th quasispecies, m^Ki = const, S is a concentration of the 
growth limiting substrate in the system, So is its concentration in the external 
flow. 

There is no apparent way to solve the system ( |27| ) and to hamiltonize it 
under arbitrary values of parameters. 

For the further analysis of this system let us take into account some qual- 
itative peculiarities of its dynamics which are known from experimental data 
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and numerical simulations [|I6fl, [|T7 



The system dynamics in general case is characterized by two basic sta- 
dium: 

the first one is an initial transient with a period of order 2/D 4/D where 
the function S(t) varies strongly and nonlinear ly; 

the second one is a quasistationary behavior in the interval 4/ D < t < r 
where r is an observation time of the system. S(t) can be approximated in 
this case by a linear function in this range up to 3-7%. Let 

S(t) =S + b(t--r). 

Here S is an average value of the S(t) during the observation time r of the 
system, b is an average velocity of the concentration variation of the rate- 
determining substrate in the system. 

Moreover, fii(S) can be also approximated by a linear function in the range 
of the quasistationary behavior where S(t) < K{. 

Expanding Monod function (|27|) in a neighborhood of the point S = S, 



we have in the second order: 

/i t (S) = a t + faS, 

A = dfn(S)/dS = rmKi/iKi + S) 2 . 
With the approximations above, the DSSF dynamics in the quasistation- 
ary stadium is described by the following equations: 

Vi = V%{cn J r fii't), (28) 
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where ol{ = an +fli(S — hbr) — D, (3{ =b • fy. It is easily to give a Hamiltonian 
form for the system flZ8|) . Choosing canonical variables in a 2N— dimensional 
phase space as 

Zi = \ogyi, a t = ai + (3i-t 

and the Hamiltonian 

N l 

h(z,a) = E(o a ? - ZiPi), 
i=l z 



we arrive to the Hamiltonian form of Eqs. 

Zi = dh/ddi = di. 



(29) 



a t = -dh/dzi = Pi. 
Note, that the hamiltonian form of the DSSF fl29| ) is similar to the DSSO 
case (|25Q and it describes uniformly accelerated rectilinear motion in space of 
variables (z^ in a constant and homogeneous external field with the potential 
U(z) = — Y,f=if3jZj determining a constant "force" F = f3. The Lagrange 
function is: 



N N I 

3=1 3=1 1 



(30) 



The variational principle has the form (|TT|) where l(z, z) is of the form ( |30| 
Canonical transformation (zi,ai) — >• (yi,ri), 

Zi = logyj, 

"i yi'ii 

allows to rewrite ( |2B| ) in Hamiltonian form, 



(31) 
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Vi = dh/dn = y?r i: 

r l = -dh/dyi = -y l r 2 + fay^ 1 . 



(32) 



with the Hamiltonian 



N l N I 

Hy,r) = Efcs/M - A logy*) = E ~g ij (y)nrj + U{y), 

i=l £ i,j'=l ^ 



^0/) = ~~ £i=iAlog|/i. The solution of the system (|52|) with the initial 
conditions y l {t)\ t =Q = yw, r t (t)\ t=0 = r i0 takes the form : 



Vi(t) = y t oexp(y l0 r l0 t + |^t 2 ), 



nit) = 



UiOTiQ + (3{t 



(33) 



2/io exp(^ no^ + \(3it 2 )' 
Expressions of yo(t, y, r), ro(£, y, r) as in the case of ([To]) , ([JJJ) are obtained 

from ( p3|) by substitution: y^ <-> y^ tq{ «-> r*, t — > — £. 

Redenoting yi ^ y\ we can write down contravariant and covariant com- 
ponents of the metrics as follows: g l ^(y) = 5^{y 1 ) 2 , gij(y) = 5ij(y l )~ 2 , respec- 
tively. The Lagrange function, / = E^Li y k Tk — h, and the Euler-Lagrange 
equations (|2T|) take here the form: 



I = \ E 9ij(yW ~ U(y), 



d 2 v J N 

• + 



+ E rLGOi/V = - E 9 ik (v) dU[v) 



N 

dt2 . . E, v., %t • 

respectively, where T J kl (y) are the Christoffel symbols related to the metrics 
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IV. DISTRIBUTION FUNCTION OF INITIAL DATA FOR DS 



Hamiltonian form of the DS dynamics allows one to investigate evolu- 
tion of an initial data distribution. Consider this problem for the DSSO 
case in more detail. Let initial data for the DSSO are given in a domain 
O of phase space A4 with a distribution function fo(p,r), fo(p,r) > 0, 
So fo{p,r)dpdr = 1. The problem of initial data indeterminancy arises for 
the Darwin systems due to natural restrictions on the exactness of measure- 
ments and Poisson character of initial values distribution of separate geno- 
types quantity. Consider the evolution of the distribution function f(t,p,r) 
under the condition that 

f(t,p,r)\ t=0 = f (p,r). (34) 

Supposing that the randomness is brought in the DS behavior only in the 
form of random distribution of initial data and further DS evolution is de- 
termined and is subjected to the Hamiltonian equations above, we deduce 
that the distribution remains constant on the phase trajectories on account 
of conservative character of the Hamiltonian dynamics. In other words, the 
distribution function f(t,p,r) obeys the Liouville equation: 

% + {f,h} = 0. (35) 



The solution of fl35| ) with the initial condition fl33Q is obtained by substitution 



of the expressions Po(t,p,r), ro(t,p, r) inverting formulas (|I5|), (|T6|) into the 
initial function _/q : 



f(t,p,r) = fo(po(t,p,r),r (t,p,r)). (36) 

As an illustration let us consider more specifically the simplest case of phase 
space N = 2 when a point of phase space A4 is defined by one coordinate 
(p) and one momentum (r) variables. In this case Po{t,p,r), r (t,p, r) are 
written down in the form: 

p (t,p,r) = ^[l-tMih(t-0 + 5)}, 

r (t,p,r) = 89 cosh 2 (t ■ 9 + 5), 
1-p) 



where 9 = \rp{l — p), 5 = \ log 



P 



f(t,p, r) = f (hl - tanh(t • 9 + 5)},89 cosh 2 (t ■ 9 + 5)). (37) 
Eq. ( |T7| ) allows one to obtain time-dependent distribution function of p: 

f(t>P) = J f{t,P,r)dr. 

Evolution of the distribution function can be inferred from the transformation 
of equiprobabilistic surfaces which are obtained from equation f(t,p,r) = 
const in different moments of time t. 

Similar results can be easily obtained in the case of DSSF. Evolution of 
the distribution function is given by ( |36| ) in which it is necessary to substitute 
the equations yo(t,y,r), ro(t,y,r) inverting (|33p . 



V. CONCLUSION 



Hamiltonian form of the dynamic equations describing Darwin selection 
process, being demonstrated here by the simple examples of the DSs, leads 
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to an interesting, by our opinion, aspects of the dynamics. For example, DS 
can be investigated in the framework of near-integrable Hamiltonian systems 
with more complicated functions D then studied in the present work. 

It is pertinent to note that in the present work the phase space M. includes 
original information space A. Such way of Hamiltonization seems to be more 
preferable if compare to one proposed in [|T8| for the Volterra-Lotka type 
systems where phase space is built up from the original variables. In the 
latter case a number of constraints are imposed on the system that restricts 
the generality. 
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